%****************************************************************
%** This program is specifically written to replicate Xtreg in **
%** Hansford and Gomez's (2010) APSR article.                  **
%** Created by Min-Hua Huang, Jan 25, 2012                     **
%** Model: Maximum Likelihhod by Fixing at Minimum, with       **
%**        Constrained Optimization, but No Correction on      **
%**        Location Estimates.                                 **
%****************************************************************
clear; clc; format shortG;
load HansfordGomez_Data.mat;
code=HansfordGomez_Data(:,2);
y=HansfordGomez_Data(:,27);
n=size(y,1); 
x0=[HansfordGomez_Data(:,26) HansfordGomez_Data(:,3) HansfordGomez_Data(:,25) HansfordGomez_Data(:,29) HansfordGomez_Data(:,11:23)];
grandmeanx0=mean(x0);
m=size(x0,2);
obs=(1:n)';
replicate=[y x0 code obs];
replicate=sortrows(replicate,19);
dmean=dmean1(replicate,19);
y2=[dmean(:,1) replicate(:,19:20)];
mean4a=[(replicate(:,1:18)-dmean) replicate(:,19:20)];
mean4a=sortrows(mean4a,20);
y=mean4a(:,1); y2=sortrows(y2,3);
y2=y2(:,1); % district mean
wfmin=ones(n,2); wf1min=ones(n,2); 
for i=1:1:2
    wfmin(:,i)=mean4a(:,i+1)-min(mean4a(:,i+1)); 
end
cr1=HansfordGomez_Data(:,10).*wfmin(:,2);
cr2=wfmin(:,1).*wfmin(:,2);
afixx=[ones(n,1) wfmin cr1 cr2 HansfordGomez_Data(:,11:23)];
ols_afix=(afixx'*afixx)\(afixx'*y);
%**** Set up initial values for constrained optimization ****
sigma=sqrt((y-afixx*ols_afix)'*(y-afixx*ols_afix)/(n-m-1));
x=afixx; 
r0=[ols_afix;sigma];
runmax=31;
dmy=13;
a=0;
b=100;
tolv=10^-4;
y1=y;
frac=1;
tau=50;
%**** Executing maximim likelihood constrained optimization ****
%*****************  Without bias correcttion *******************
fs=ML1(y1,x,r0,a,b,runmax,dmy,tolv,tau,y2,frac); 
ml1b=fs(1:m+1)';
ml1sigma=fs(m+2);
ml1se=robuster(y1,afixx,[ml1b;ml1sigma],code);
ml1t=ml1b./ml1se;
hr=boundary(y2,frac); 
MLCO1ymax=max(x*ml1b+y2);
MLCO1ymin=min(x*ml1b+y2);
xmax=zeros(m+1,1); xmin=zeros(m+1,1);
   for i=1:1:m+1
       xmax(i)=max(x(:,i));
       xmin(i)=min(x(:,i));
   end    
xlimit=[xmax xmin];
v=vsign(ml1b);
MLCO1ci=cpar(fs(1:m+2)',a-hr(2),b-hr(1),x,xlimit,v,dmy);
MLCO1bv=sum(MLCO1ci>tolv);
MLCO1predict=range1(ml1b,hr(2),hr(1),x,xlimit,v,dmy,y2);
fprintf('\nReplication of the column 1 in Table 1, Hansford and Gomez (2010)\n')
fprintf('\nMaximum Likelihod Estimation by Fixing at Minimum with Constrained Optimization\n')
fprintf('This Model Does Not Have Bias Correction in the Demeaning Operation\n')
fprintf('                                 beta      se       t             \n');
fprintf('Partisan Composition           % 7.3f  % 7.3f  % 7.3f\n',ml1b(2),ml1se(2),ml1t(2));
fprintf('Turnout                        % 7.3f  % 7.3f  % 7.3f\n',ml1b(3),ml1se(3),ml1t(3));
fprintf('Turnout*GOP Incumbent          % 7.3f  % 7.3f  % 7.3f\n',ml1b(4),ml1se(4),ml1t(4));
fprintf('Turnout*Partisan Composition   % 7.3f  % 7.3f  % 7.3f\n',ml1b(5),ml1se(5),ml1t(5));
fprintf('Yr52                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(6),ml1se(6),ml1t(6));
fprintf('Yr56                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(7),ml1se(7),ml1t(7));
fprintf('Yr60                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(8),ml1se(8),ml1t(8));
fprintf('Yr64                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(9),ml1se(9),ml1t(9));
fprintf('Yr68                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(10),ml1se(10),ml1t(10));
fprintf('Yr72                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(11),ml1se(11),ml1t(11));
fprintf('Yr76                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(12),ml1se(12),ml1t(12));
fprintf('Yr80                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(13),ml1se(13),ml1t(13));
fprintf('Yr84                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(14),ml1se(14),ml1t(14));
fprintf('Yr88                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(15),ml1se(15),ml1t(15));
fprintf('Yr92                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(16),ml1se(16),ml1t(16));
fprintf('Yr96                           % 7.3f  % 7.3f  % 7.3f\n',ml1b(17),ml1se(17),ml1t(17));
fprintf('Yr2000                         % 7.3f  % 7.3f  % 7.3f\n',ml1b(18),ml1se(18),ml1t(18));
fprintf('Constant                       % 7.3f  % 7.3f  % 7.3f\n',ml1b(1),ml1se(1),ml1t(1));
fprintf('Sigma                          % 7.3f\n',ml1sigma);
fprintf('Max yhat                       % 7.3f\n',MLCO1ymax);
fprintf('Min yhat                       % 7.3f\n',MLCO1ymin);
fprintf('Theoretical Boundary Violations% 7.0f\n', MLCO1bv);
fprintf('N                              % 7.0f\n', n);
fprintf('\nPredicted Value by Varying        Min      Max\n');
fprintf('All Variables                  % 7.1f  % 7.1f\n',MLCO1predict(1,1),MLCO1predict(1,2));
fprintf('Constant                       % 7.1f  % 7.1f\n',MLCO1predict(2,1),MLCO1predict(2,2));
fprintf('Partisan Composition           % 7.1f  % 7.1f\n',MLCO1predict(3,1),MLCO1predict(3,2));
fprintf('Turnout                        % 7.1f  % 7.1f\n',MLCO1predict(4,1),MLCO1predict(4,2));
fprintf('Turnout*GOP Incumbent          % 7.1f  % 7.1f\n',MLCO1predict(5,1),MLCO1predict(5,2));
fprintf('Turnout*Partisan Composition   % 7.1f  % 7.1f\n',MLCO1predict(6,1),MLCO1predict(6,2));
fprintf('Yr52                           % 7.1f  % 7.1f\n',MLCO1predict(7,1),MLCO1predict(7,2));
fprintf('Yr56                           % 7.1f  % 7.1f\n',MLCO1predict(8,1),MLCO1predict(8,2));
fprintf('Yr60                           % 7.1f  % 7.1f\n',MLCO1predict(9,1),MLCO1predict(9,2));
fprintf('Yr64                           % 7.1f  % 7.1f\n',MLCO1predict(10,1),MLCO1predict(10,2));
fprintf('Yr68                           % 7.1f  % 7.1f\n',MLCO1predict(11,1),MLCO1predict(11,2));
fprintf('Yr72                           % 7.1f  % 7.1f\n',MLCO1predict(12,1),MLCO1predict(12,2));
fprintf('Yr76                           % 7.1f  % 7.1f\n',MLCO1predict(13,1),MLCO1predict(13,2));
fprintf('Yr80                           % 7.1f  % 7.1f\n',MLCO1predict(14,1),MLCO1predict(14,2));
fprintf('Yr84                           % 7.1f  % 7.1f\n',MLCO1predict(15,1),MLCO1predict(15,2));
fprintf('Yr88                           % 7.1f  % 7.1f\n',MLCO1predict(16,1),MLCO1predict(16,2));
fprintf('Yr92                           % 7.1f  % 7.1f\n',MLCO1predict(17,1),MLCO1predict(17,2));
fprintf('Yr96                           % 7.1f  % 7.1f\n',MLCO1predict(18,1),MLCO1predict(18,2));
fprintf('Yr2000                         % 7.1f  % 7.1f\n\n',MLCO1predict(19,1),MLCO1predict(19,2));
